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Abstract: The local density of states and response to an incident plane 
wave of a finite sized photonic crystal (PC) with nonlinear material (NLM) 
is analyzed. Of particular interest is the excitation of surface wave modes at 
the truncated surface of the PC, which is collocated with the NLM material. 
We compute the 2D Green function of the PC with linear material and then 
include the Kerr NLM in a self-consistent manner. The 2D PC consists of 
a square array of circular rods where one row of the rods is semi-circular 
in order to move the surface wave defect mode frequency into the band 
gap. Since the surface modes are resonant at the interface, the NLM should 
experience at least an order of magnitude increase in field intensity. This 
is a possible means of increasing the efficiency of the PC as a frequency 
conversion device. 
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1. Introduction 

The refinement of Hthographic and other fabrication techniques used in photonic crystal (PC) 
development has generated a broad range of experimental, and theoretical interest [1, 2, 3, 4, 
5, 6, 7]. In particular, the ability of PC structures to control the propagation of electromagnetic 
fields has stimulated certain emerging concepts and devices that are based on PC technology, 
including microscopic lasers, optical switches, light localization, and resonant cavity antennas 
[4]. 

It is computationally convenient when investigating PC structures, to assume the PC is of 
infinite extent. This means a perfect crystal, or a crystal with defects in a supercell configura- 
tion, and this works well for band structure calculations. This approximation has limitations, 
since experimental situations involving finite-sized crystals may involve defect modes that re- 
sult from truncating the infinite periodicity, or reflections and interference effects of the fields 
around the boundaries [2]. The vast majority of theoretical studies also assume the PC con- 
sists of linear material. When characterizing the electromagnetic properties of PC structures, 
the transmittance characteristics often yield limited information regarding the band gap, owing 
mainly to the directional dependence of the incident field [6] . The quantity of experimental rel- 
evance, that affords a more complete and fundamental picture of the spectra and modes of the 
system, is the spatially and energy resolved local density of states (LD0S).[7] 

In this paper, we compute LDOS and field characteristics for a finite two-dimensional PC 
consisting of both linear and nonlinear material. Some researchers [8] have focused on surface 
modes associated with a one-dimensional PC. Our emphasis is on the properties of surface res- 
onant modes as a localized field to concentrate the electric field in the presence of NLM. For 
the case of two-dimensional truncated structures consisting of an array of cylinders, surface 
electromagnetic waves can be supported. These surface modes are localized near the trunca- 
tion boundary and propagate parallel to the interface. Since the intensity of the surface modes 
is concentrated at the interface, and if the interface contains nonlinear material (NLM), then 
a more efficient frequency conversion might be achieved. This could be useful in frequency 
conversion devices for use in detectors to operate between visible and infrared wavelengths. 

In Sec. 2, we outline the integral formulation of the problem. Next, we briefly discuss the 
iterative numerical algorithm used to calculate the electromagnetic fields and Green's function, 
which gives essential information regarding the LDOS. The LDOS will help reveal certain 
defect modes of the PC that occur at frequencies within a band gap. We present the results in 
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Fig. 1. Schematic of photonic crystal (PC) cross-section. The PC consists of infinitely long 
rods in a 2D array of 7 rods x 6 rods. Note that the rods in the topmost layer have been 
truncated along their diameter into a semi-circular shape. This truncation will allow local- 
ized surface modes to exist along this side of the PC within band gap frequencies. The 
lines AC and BD are guides for the eye in relating to numerical results below. The arrow 
and angle 6 denote an incident plane wave. The PC has period a, the cylinders have radius 
r = 0.2a with permittivity £c = 8.9, and the background medium is vacuum. Dispersion 
characteristics of this PC are given in Fig. 2. 



Sec. 3 and find that within the PC structure, the LDOS is small but nonzero within the band gap, 
and that any localized states reside mainly within the truncated surface. Results are also shown 
that illustrate the effect of incorporating NLM material. 

2. Theory 

We calculate the LDOS of a finite size 2D PC for the cases where the dielectric material is all 
linear or some rods are nonlinear. We will also investigate the response of the PC to an incident 
electric field. Of particular interest are localized surface wave modes along a PC interface. To 
model the 2D PC, we consider a finite- sized square array of dielectric cylinders as shown in 
Fig. 1. We begin with the case where the PC consists of linear material and assume the electric 
vector is polarized parallel to the rods. Since the PC considered here is finite size, the possibility 
of defect and localized modes exist. The design parameters of the PC are such that a band gap 
exists and truncation of the topmost layer of rods permits localized surface wave modes to occur 
at frequencies within the band gap. The dispersion of the bulk PC and surface modes is shown 
in Fig. 2. 

In the case of linear media, the LDOS is independent of an incident field. In the case of 
non-linear material, we illuminate the PC with a plane wave and include non-linear material 
effects in a self-consistent manner. In this case, any change in the LDOS as calculated here, is 
examined. We neglect the effects of coupling to higher harmonic fields. 

2.1. Dyson 's equation and iterative solution 

The method of analysis involves iterative solution to Dyson's equation to calculate the Green 
function and the Lippman-Sch winger equation for the electric field. [10, 11, 12] We assume 
harmonic fields that vary as exp(— /(O^). Here we consider an especially simple case where the 
electric vector is polarized parallel to the axis of the dielectric cylinders. In this case, the Green 




Fig. 2. (a) Band diagram of the bulk PC for an infinite 2D array of round rods with the 
electric vector parallel to the rods (E||) (blue) and the magnetic vector parallel to the rods 
(H||) (red). The vertical axis is normalized frequency / = {(o/c){a/2n). As shown by 
the blue curves, an infinite PC of all round rods with the parameters given in Fig. 1 has a 
complete band gap in the frequency range 0.32 </ < 0.44 for (E| | ) modes. For (H| | ) modes, 
there is no complete bandgap. (b) These data are a subset of the F-X portion of (a) where 
the green curve shows the dispersion of E| | surface modes that propagate with wavevector k 
parallel to the semi-circular cylinder interface at frequencies within the bandgap. The black 
line is the light line and the gray lines show the band gap limits. Of particular interest in 
the analysis to follow is the green curve in (b) for the localized surface wave mode. We 
have chosen to truncate the top row of cylinders along their diameter, however varying the 
amount of truncation could be done and this would yield different dispersion curves for the 
surface mode. The data in (a) and (b) are reproduced from [9]. 

dyadic G = zzG^z effectively becomes a scalar that satisfies Dyson's equation 

G,,{p,p') =gUp,P')+i4 ld^p''Gi{p,p'')Hp'')G,,{p'',p') (1) 

and the electric field E = zE^ satisfies either of the following forms of the Lippmann-Schwinger 
equation, 

E,{p)=E'^{p)+kl [ d^p'G^{p,p')e{p')Et{p') (2) 
E,{p)=E'^{p)+kl [ d^p'Gi{p,p')HP')Ez{p') (3) 

J A 

where ko = co/c and p = (x^y). The integration in Eqs. (l)-(3) is only over the cross-sectional 
area of the cylinders. The G^^{p^p') = ^Hq {ko\p — p'\) is the Green function for free space 
in the absence of the PC structure and G^z is the Green function with the effects of the PC 
included. The E^{p) is the incident electric field. Note that when p = p' in Eqs. (l)-(3) the 
singularity associated with the Green function would normally require special treatment, but in 
the simple case treated here, the singularity is more easily handled. [13] For simplicity, we omit 
all z subscripts henceforth. 

The PC structure is described by the £ (p) that is non-zero only within the dielectric rods. 
Initially, we assume the rods are homogeneous with permittivity £c embedded in a background 
medium of unit permittivity. This yields 



e(p) = ec-l 



(4) 



when p is within any rod and zero otherwise. Numerical solution to Eqs. (1) and (3) can be 
performed by direct inversion, but this would typically be very large in terms of memory re- 
quirements, especially in Eq. (1). Here we provide only a brief description of the iterative 
method. The computational space is divided into an array of grids where the ith grid has area 
AA. The center of the ith and jth grid is located at pi and pj. The space is initially empty and the 
PC is then constructed by adding one "infinitesimal" £ element at a time. When the nth element 
is added we have 

G'}j = Glr^+AGl-'et.Glj i^k„, j^K (5) 
gF.„ = (/-K4„)"'g«-i i^K, j = K 

where /:qAA = A. To initiate the iteration algorithm, set n = 1 and there is a single perturbation 
Ej^^ at grid k\ . The term represents the integration term when / = j = k\ as 



< = / gO , dA = - 1 (6) 



inkoRHi{koR) 



where R = ^J bAjn. In (1), the integral over the grid square has been approximated by integrat- 
ing over an equivalent circular area of radius R. The electric field solution can be built up in a 
similar iterative fashion or calculated as 

£, =£f+Af; Gi^„4„< (7) 

where is the incident field at the /th coordinate and N is the number of perturbation elements 
E]^^ in the structure. We have given a very brief outline of the method since details[10, 11, 12] 
are described elsewhere and the method is not the main focus of this work. 

2.2. Local density of states 

Computation of the electric field and the Green dyadic yields much information, including the 
normalized LDOS, that is defined here as 

- 3[G0(p,p)] ' 

where is the free space Green function as given in Eq. (1), and following Eq. (3), and 3 
indicates the imaginary part. In the case of linear media, the LDOS as calculated from Eq. (8) 
is only dependent on the material and structural parameters of the PC. If nonlinear material 
is introduced, then obviously an excitation field is required to activate the nonlinearity. In this 
case, we still use Eq. (8) to calculate the LDOS with the PC illuminated by an incident plane 
wave of unit amplitude and at a given angle of incidence. This is done self-consistently yielding 
changes in the permittivity as described in the next section. The resulting changes in permit- 
tivity yields different values for G(p,p) and consequently the LDOS(p). The modeling of the 
nonlinearity is discussed in the following section. 

In [7], the local density of states is written as LDOS = — ^3 [G(p , p)] . Since in this work we 
consider the LDOS in Eq. (8) that is normalized to free space, we note that, — 3G^(p , p) = 0.25. 
It follows that the normalized LDOS given in numerical results below can easily be compared 
to free space relative to 0.25. 
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Fig. 3. (a) Normalized LDOS(x) vs position x within the PC for various frequencies with 
Unear material and (b) Normalized LDOS(x) at resonance frequency compared with non- 
linear material. In (a), and referring to Fig. 1, the LDOS(x) is shown on a linear scale and is 
calculated along line AC (between rods) and BD (through rods) where the BD data has been 
displaced one unit upward for clarity. For reference, the large black circles in the lower part 
of the figures indicate the positions of the rods in relation to the LDOS. In (a), the dimen- 
sionless frequency parameter / = {(o/c){a/27i) is varied for (O within the band gap from 
0.34 to 0.36 with a peak at / = 0.344. Along line BD, there is clear indication of local- 
ized states concentrated at the semi-circular rod interface {x ^ 130A, where A is the spatial 
resolution) whereas no localization at the opposite circular rod interface (x ^ lOA). Along 
line AC there is no comparable indication of localized states. It follows that the localized 
states exclusively reside within the semi-circular rods. This indicates the importance of de- 
fect parameters, such as altering the rod cross-section, in the surface mode dispersion. In 
(b), the resonant LDOS(x) is plotted on a log scale along lines AC and BD for / = 0.344 
with X = ^ ^rid 0.01. The X = ^ curves are repeated from (a). To activate the NLM, the PC 
is illuminated by a plane wave of unit amplitude incident at = 45^. Along line AC, the 
green and black curves are superimposed. Along line BD, the red and blue curves are nearly 
superimposed with the blue curve slightly larger in value that indicates a slight reduction in 
LDOS with non-zero X- This reduction is not very clear on this log scale, however the log 
scale does show more clearly the LDOS(x) throughout the PC. This indicates that the pres- 
ence of non-linear material in the semi-circular rods does not significantly alter the LDOS. 
Finally, these numerical data are completely independent of an incident field. 



2.3. Non- linear material 

Others have considered the effect of NLM on the infinite photonic crystal. [14, 15, 16, 17, 18] 
Using the iterative procedure[10, 11, 12] briefly described above, we calculate for a finite size 
PC the Green function, LDOS, and the response of the PC shown in Fig. 1 to an incident plane 
wave. To consider the addition of a NLM in the topmost layer of semi-circular rods, we take a 
self-consistent approach. [15] Assuming a Kerr non-linearity we write 

{Sc — I -\- x\Ec{p)\^ '• P within any semicircular cylinder 
£c — 1 : p within any round cylinder (9) 
: p otherwise 

Initially, the electric field distribution within the semi-circular rods, E^{p), is obtained when 
£^{p) = e^-l. Using E^{p), the permittivity is then changed to e^(p) = - 1 +;^|^c (P)P- 
Using (p), a new Green function Gij and field distribution E}. (p) are calculated. This process 



Plane Wave Incident on PC at /= 0.344 and = 45 
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Plane Wave incident on PC at / = 0.344 for 6=135 " 
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Fig. 4. Plane wave at normalized frequency / = 0.344 incident on PC at angles of (a) 45^ 
and (b) 135^ (see Fig. 1) for both linear and non-linear material. The x-coordinate is in 
A units and the electric field is plotted on a log scale as \E\^ since this is relevant to a 
Kerr nonlinearity. For both angles of incidence, there is significant excitation of the surface 
defect mode as noted by the peaks around x ^ 130A and lines BD . In (a), there are incident 
and reflected fields in the x < PC interface region. In (b), there are again the incident and 
reflected fields for x > 130A in addition to the surface mode excitation field. Note also in 
(b) that the low field intensity for x region is because this is a near-field shadow region. 
It is seen that for a unit amplitude incident wave, the localized \E\^ field is well in excess 
of an order of magnitude larger because of the surface mode excitation and the localization 
is precisely at the sites of the nonlinear material. 
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Fig. 5. Contour plots of a plane wave at frequency / = 0.344 incident on PC at angles 
of (a) 45^ and (b) 135^ for linear material. The angles of incidence are indicated by the 
arrows. These plots show the electric field intensity throughout the PC and within 2 periods 
a outside the PC, The line plots of Fig. 4 are a subset of these data. The x and ^-coordinates 
are in Xq = a/0.35 units. Surface defect mode excitation is clearly seen along the truncated 
layer of rods. Note that the outline of the PC rod structure is superimposed on this contour 
plot. 



is repeated until succeeding electric fields throughout the PC converge to a stable result. This 
yields not only the electric fields, but also any changes in the LDOS. 

In testing for convergence of the self-consistent algorithm, we found that for % < 0.01 ap- 
proximately 10 iterations or less were needed to converge. If the discretized computational 
domain contains M points we calculate the electric field E^^\x,y) for all M points. 

We assume convergence is satisfied when j^T^xj I ^^^\^^y) |^ 0.0001 where A^\x^y) = 
[E^^\x,y) — E^'^~^\x^y)]/ E''^~^\x^y) and n is the number of iterations. We did not do an ex- 
tensive parametric evaluation of the self-consistent convergence conditions. We did note that 
when ;^ = 0.1, the self-consistent algorithm failed to converge. In any case, we did not try to 
use experimentally obtained X values, rather we note that X ^ 0-01 is well within any realistic 
material values. 

3. Numerical Results 

The PC is modeled as a square array of dielectric rods with period a and radius r with physical 
parameters [9] chosen so that there is a complete band gap for the infinite crystal and this is 
shown in Fig. 2(a) and the defect mode surface wave dispersion is shown in Fig. 2(b). The 
computational space is divided into square grids of side dimension A and the PC period is 
a = NA. We choose a frequency cOq within the band gap where {(O^/ c){a/27t) = 0.350. We 
set N = 20 and this yields a free space resolution of A = 0.0175Ao with cOq/c = 271 /Xq. The 
resolution in a cylinder is then degraded to A^ = AVK9, and this is still adequate for accurate 
computation. 

To check the accuracy of our results, we have compared calculation of the electric field as 
obtained from Eqs. (2) and (3). We note that for reasonable sized grid of N points, Eq. (2) can 
be solved for Ez{p) by direct inversion. Also, Ez{p) can be obtained from Eq. (3) after first 
obtaining G^^{p^p'). We find that the electric fields computed by both methods are in excellent 
agreement and we conclude that the method of calculating G^(p,p') is also accurate. 

Figure 3 shows the results of computations of the LDOS associated with the PC and the plane 
wave response of the PC. In Fig. 3(a), the calculated LDOS is shown for several frequencies 
and a PC with all linear material. At the surface mode resonant frequency, and with the PC il- 
luminated with a unit- amplitude plane wave at = 45^, the computed effect of including NLM 
is shown in Fig. 3(b). The NLM coefficient is chosen as ;^ = 0.01. The results of the LDOS cal- 
culations indicate a large concentration of states available at the truncated-rod interface. When 
the NLM coefficient x is non-zero, there is noticeable decrease in the LDOS. This decrease is 
due to the change in the permittivity of the semi-circular rods and the surface wave mode is 
quite sensitive to any physical or material changes in the semi-circular interface. 

In Fig. 4 we consider the response of the PC when illuminated by a plane wave at two 
different angles relative to the truncated rod interface. We plot the electric field as \E{x^y)\^ 
since this is the magnitude relevant to Kerr nonlinearity. For both angles of incidence, there 
is strong excitation of the surface mode. It is seen that for ;^ =0.0 and 0.01 for both angles 
of incidence B = 45^ and 135^, the surface mode peaks are quite similar in both cases. This is 
likely due to the fact that very little of the PC actually contains nonlinear material. Nevertheless, 
there are some differences in the \E{x^y)\^ results when x = 0-0 and 0.01. The addition of 
;^ = 0.01 causes a slight decrease in \E{x^y)\^ and this is consistent with the LDOS data in Fig. 
4(b). 

A more comprehensive view of the field profile is illustrated in Fig. 5, which shows the 
corresponding contour plots relating to Fig. 4 for linear media. Here the field magnitude is 
plotted over the entire PC including two periods outside the PC boundary. The surface mode 
excitation at the semicircular rod sites is clearly seen. The field magnitude inside the PC is quite 
diminished and the location of the rods comprising the PC is shown by an overlay of an outline 



of the rods. 



4. Conclusions 

An analysis of the LDOS for surface wave modes on a truncated 2D photonic crystal interface 
has been accomplished using a Green function formalism. The PC is composed of circular rods 
(except for the semicircular interface) and the electric field is parallel to the rods. The results 
show that the presence of NLM does not significantly alter the LDOS as defined in this paper 
and this could be due to the NLM being a small fraction of the total volume of PC material. 
If the idea of using a 2D PC as a component in frequency conversion is viable, then it seems 
that additional volume of NLM would probably be needed. We have also neglected the fact that 
frequency conversion does occur and that there is an additional coupled equation to consider 
that would describe the generation of higher harmonic fields due to the Kerr nonlinearity. This 
may be a source of damping that could broaden the resonance associated with the surface wave 
mode. This could reduce the conversion efficiency. 
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